Skip to content

Regenerate the pH 7 protonation bundles with Marvin 26.1 - #295

Merged
samseaver merged 1 commit into
ModelSEED:devfrom
freiburgermsu:marvin-26.1-pka-regeneration
Sep 21, 2026
Merged

samseaver merged 1 commit into
ModelSEED:devfrom
freiburgermsu:marvin-26.1-pka-regeneration

Conversation

@freiburgermsu

Copy link
Copy Markdown
Member

Closes the item Marvin_26.1_pKa_Run.md left under "Not done here": protonations/ still held only the 23.4 bundle. Companion to #289 / #292.

53,127 compounds, 141,715 rows — the row count 23.4 shipped, exactly, in every source. ChEBI, KEGG and Rhea cover an identical compound set. MetaCyc swaps one in each direction. Zero unparseable, zero non-standard InChI, zero InChIKeys disagreeing with the InChI beside them.

source compounds rows 23.4 rows
ChEBI 11,544 30,488 30,488
KEGG 16,275 46,911 46,911
MetaCyc 25,071 63,665 63,665
Rhea 237 651 651
total 53,127 141,715 141,715

Two things to weigh before merging

1. The tautomer step is missing

ChargeMol.java took the dominant tautomer at pH 7 before the major microspecies. That needs TautomerizationPlugin, licensed under the Isomers Plugin Group; these regenerations run on a Protonation-group licence. Verified, not assumed:

$ cxcalc "CC(=O)CC(=O)C" majortautomer -H 7
chemaxon.license.api.LicenseException: No valid license has been found.
Product name: Isomers Plugin Group

Net charge at pH 7 is unchanged for 79.7% of shared compounds and different for 20.3%, skewed toward deprotonation (11.8% one unit lower, 4.1% two or more; 3.6% one higher).

That 20.3% is an upper bound on the tautomer effect, not a measurement of it. It also contains the genuine 23.4 → 26.1 engine change, which the pKa work already showed moves site counts on roughly one shared set in seven — a pKa crossing 7 flips the protonation state on its own. The two cannot be separated without an Isomers licence, and the report does not pretend otherwise. Restoring the step is a few lines if that licence is obtained.

2. Marvin's InChI and SMILES importers lose stereochemistry

Likely the most reusable finding here, and not documented anywhere obvious. Round-tripping the source InChI with no protonation at all — import, export, compare:

import path reproduces the source InChI
MolImporter.importMol(<InChI>) 61%
MolImporter.importMol(<SMILES>) 39%
MolImporter.importMol(<RDKit molblock>) 99.6%
RDKit alone (control) 99.7%

The loss is on import, not export — what goes is the /t and /m parity layers. Every structure is therefore bridged through RDKit into a molblock before Marvin sees it, which is what ChargeMol.java was doing when it read args[0] as a mol file rather than a string.

Measured on the compounds where protonation changes nothing, so output should equal input:

source no-op compounds this run 23.4
ChEBI 7,128 99.0% 96.6%
KEGG 8,346 99.0% 97.8%
MetaCyc 13,738 98.6% 96.0%
Rhea 135 98.5% 98.5%

Unbridged this scores 55–61% — and every coverage count still looks perfect. inchi:AuxNone,SAbs on export matters just as much: without SAbs, 30.9%.

Smaller notes

  • InChIKey is derived from the InChI actually written. Marvin's inchikey export ignores SAbs and disagrees with its own inchi export — for POLYMER_10033 it emits ...-ARPYZQPTNA-N where the InChI it just wrote hashes to ...-VFUOTHLCSA-N, which is what 23.4 shipped. 0 disagreements across all 44,294 InChI rows.
  • One engine, not the CLI. Plugin and cxcalc majorms agree on 296/300 ChEBI compounds; all 4 exceptions are organometallics where cxcalc silently drops the metal fragment (C35H34MgN4O5C35H34N4O5). Same reason ChargeMol.java fragmented and re-fused.
  • Query molecules are recovered, not dropped. aromatize() on query atoms leaves bonds no SMILES writer can express, and MetaCyc spells its protein-bound cofactors with ~ any-bonds. Without handling both, 156 MetaCyc compounds would have been lost silently.
  • Charge on wildcard-adjacent atoms is inherited. 23.4 wrote *[NH2+][C@@H](... for POLYMER_12621; so does this. Correcting it is a curation decision about what * means, so it is recorded and left alone.
  • Failures, by id. Four compounds yield no row; three are H⁺ (InChI=1S/p+1, no heavy atom — 23.4 has no row either, so agreeing is correct). The one real failure is MetaCyc HypC-Dimer-Fe-CO2, an ArrayIndexOutOfBoundsException in PkaPlugin.run() — one row out of 141,715.
  • Re-running rewrites a bundle byte for byte.

Not done here

  • Compound records untouched — Update_Compound_Structures_Formulas_Charge.py is the step that rewrites them.
  • Print_Structure_Formula_Charge.py has not been re-run over these files; the formula/charge columns are Marvin's own.
  • ChEBI keeps bare ids in this bundle, matching the 23.4 protonation file (the pKa bundles use CHEBI_). The id-format migration stays separate.

Full detail in Biochemistry/Structures/_reports/Marvin_26.1_Protonation_Run.md.


Generated by Claude Code

Closes the item Marvin_26.1_pKa_Run.md left under "Not done here":
protonations/ still held only the 23.4 bundle.

53,127 compounds, 141,715 rows -- the row count 23.4 shipped, exactly, in
every source. ChEBI, KEGG and Rhea cover an identical compound set; MetaCyc
swaps one (loses HypC-Dimer-Fe-CO2 to a PkaPlugin crash, gains
A-DNA-WITH-OPPOSING-AP-SITE). Zero unparseable, zero non-standard InChI,
zero InChIKeys disagreeing with the InChI beside them.

Two things a reviewer should weigh before merging.

THE TAUTOMER STEP IS MISSING. ChargeMol.java took the dominant tautomer at
pH 7 before the major microspecies. That needs TautomerizationPlugin, which
is licensed under the Isomers Plugin Group; these regenerations run on a
Protonation-group licence, and cxcalc majortautomer returns a
LicenseException. Net charge is unchanged for 79.7% of shared compounds and
different for 20.3%, skewed toward deprotonation. That 20.3% is an upper
bound on the tautomer effect, not a measurement of it -- it also contains
the genuine 23.4 -> 26.1 engine change, which the pKa work already showed
moves site counts on one shared set in seven. The two cannot be separated
without that licence. Restoring the step is a few lines if it is obtained.

MARVIN'S InChI AND SMILES IMPORTERS LOSE STEREOCHEMISTRY. Round-tripping
the source InChI with no protonation at all reproduces it 61% of the time
through the InChI importer and 39% through SMILES, against 99.6% through an
RDKit-built molblock and 99.7% for RDKit alone. The loss is on import, not
export. Every structure is therefore bridged through RDKit into a molblock
before Marvin sees it -- which is what ChargeMol.java was doing when it read
args[0] as a mol file rather than a string. On the compounds where
protonation changes nothing, and the output should equal the input, this run
scores 98.6-99.0% against 23.4's 96.0-98.5%; unbridged it scores 55-61% and
every coverage count still looks perfect.

Also: the InChIKey is derived from the InChI actually written, because
Marvin's inchikey export ignores SAbs and disagrees with its own inchi
export. Charge on wildcard-adjacent atoms is inherited from 23.4, not
introduced, and is left alone as a curation question. Re-running rewrites a
bundle byte for byte.

Compound records are untouched; Update_Compound_Structures_Formulas_Charge.py
is the step that rewrites them, and Print_Structure_Formula_Charge.py has not
been re-run over these files.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@samseaver
samseaver merged commit 41b20c2 into ModelSEED:dev Sep 21, 2026
freiburgermsu added a commit to freiburgermsu/ModelSEEDDatabase that referenced this pull request Sep 22, 2026
Review of ModelSEED#295 found compounds losing their R groups. Confirmed, and the
cause was this bundle's formula column, not Marvin and not the structures.

WHAT WAS WRONG. The first cut wrote Marvin's getFormula() into the formula
column and deferred Print_Structure_Formula_Charge.py as a follow-up step.
Marvin omits wildcard atoms from a formula; this repository renders them as R,
in one line of that script:

    formula = re.sub(r'\*', 'R', formula)

So SMILE rows whose formula contains R fell from 23.4's 8,712 to 8.
Stearoyl-ACPs went from C32H60N3O9PR2S to C32H60N3O9PS, and
Update_Compound_Structures_Formulas_Charge.py carried that into 6,052 compound
records -- cpd00049 "carboxylic acid" became CHO2, a generic compound quietly
ceasing to be generic.

THE STRUCTURES WERE NEVER AFFECTED: 8,727 SMILE structures carry a `*` in both
bundles, identically. Only the column was wrong. That is why it survived the
whole validation suite -- coverage, compound sets, InChI agreement, InChIKey
consistency and stereochemical fidelity all inspect structures, and the
structures were always right. A column-only defect was invisible to every
check. A formula-column comparison against 23.4 is now part of the report.

THE FIX. formula and charge come from Print_Structure_Formula_Charge's own
parse_structure -- imported, not reimplemented -- computed per row from that
row's structure string. Deferring it was the error; the cascade consumes this
column directly, so the bundle has to be correct as written.

Two further cases had to be reconciled to make that hold everywhere:

  parse_structure's substitution only fires when RDKit produced the formula,
  since RDKit renders a dummy atom as `*`. OpenBabel and Marvin both omit them,
  so structures like ISOCITHASE-P's `*OP(=O)(=O)=O` -- an invalid valence only
  OpenBabel will read -- still lost their R. The convention is now enforced as
  an invariant: a SMILE structure carrying `*` gets an R in its formula. This
  bundle violates it 0 times; 23.4 violated it 23 times.

  Counting wildcards to enforce that has its own trap: a molblock R atom, which
  is how RDKit writes every dummy atom and therefore how every structure
  arrives through the import bridge, reads back from Marvin as symbol "R#", not
  "R". A wildcard set without R# counts zero on a molecule that plainly has
  them.

RESULT. R-bearing SMILE rows 8,735 against 23.4's 8,712; per source ChEBI
2,071=2,071, KEGG 964=964, Rhea 30=30, MetaCyc 5,670 vs 5,647. Formula-column
agreement with 23.4 rises from 68,466 rows to 76,479. Coverage, structures,
InChI, InChIKey and stereochemistry are unchanged -- 141,715 rows, still
exactly 23.4's count, 0 non-standard InChI, 0 InChIKey disagreements.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
freiburgermsu added a commit to freiburgermsu/ModelSEEDDatabase that referenced this pull request Sep 22, 2026
Self-review of the numbers in the previous commit, not a data change --
the bundles are untouched and were already correct.

Counting R groups with `"R" in formula` also matches the elements Ru, Rb, Rh,
Re and Rn. That inflated every figure by 8 and, worse, made the regression look
milder than it was: the eight rows that appeared to survive with an R group
were ruthenium and rubidium compounds, not R groups at all.

  SMILE rows with an R group      was reported   actually
    this bundle                        8,735       8,727
    marvin 23.4                        8,712       8,704
    as merged in ModelSEED#295                      8           0

So ModelSEED#295 did not leave 8 R groups standing; it left none. Per source the
corrected counts are ChEBI 2,071 = 2,071, KEGG 957 = 957, Rhea 30 = 30, and
MetaCyc 5,669 against 5,646 -- the +23 delta, and every other conclusion, is
unchanged.

Re-verified against the committed bundles while correcting this:

  coverage         141,715 rows and 53,127 compounds, both exactly 23.4's
  invariant        0 violations in either direction; 23.4 has 23
  vs merged ModelSEED#295   9,185 formulas and 84 charges changed, 0 structures,
                   0 rows added or removed
  columns          one generated_on, tool_version and ph across all four files
  reproducibility  re-running a source rewrites it byte for byte

Both the report and the script docstring now state the rule (`R(?![a-z])`) next
to the numbers, so the next person to check does not repeat it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
freiburgermsu added a commit to freiburgermsu/ModelSEEDDatabase that referenced this pull request Sep 22, 2026
Self-review of the numbers in the previous commit, not a data change --
the bundles are untouched and were already correct.

Counting R groups with `"R" in formula` also matches the elements Ru, Rb, Rh,
Re and Rn. That inflated every figure by 8 and, worse, made the regression look
milder than it was: the eight rows that appeared to survive with an R group
were ruthenium and rubidium compounds, not R groups at all.

  SMILE rows with an R group      was reported   actually
    this bundle                        8,735       8,727
    marvin 23.4                        8,712       8,704
    as merged in ModelSEED#295                      8           0

So ModelSEED#295 did not leave 8 R groups standing; it left none. Per source the
corrected counts are ChEBI 2,071 = 2,071, KEGG 957 = 957, Rhea 30 = 30, and
MetaCyc 5,669 against 5,646 -- the +23 delta, and every other conclusion, is
unchanged.

Re-verified against the committed bundles while correcting this:

  coverage         141,715 rows and 53,127 compounds, both exactly 23.4's
  invariant        0 violations in either direction; 23.4 has 23
  vs merged ModelSEED#295   9,185 formulas and 84 charges changed, 0 structures,
                   0 rows added or removed
  columns          one generated_on, tool_version and ph across all four files
  reproducibility  re-running a source rewrites it byte for byte

Both the report and the script docstring now state the rule (`R(?![a-z])`) next
to the numbers, so the next person to check does not repeat it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
freiburgermsu added a commit to freiburgermsu/ModelSEEDDatabase that referenced this pull request Sep 22, 2026
Review of ModelSEED#296 flagged 74 formula rows regressing on metals. Chasing it found
the formula was a symptom. The structures were broken.

    source  Cl[Sn](c1ccccc1)(c1ccccc1)c1ccccc1      intact
    23.4    Cl[Sn](c1ccccc1)(c1ccccc1)c1ccccc1      correct
    ModelSEED#295    [Cl-].[SnH3+].[c]1ccccc1.[c]1ccccc1.[c]1ccccc1

INCHI DISCONNECTS METAL-LIGAND BONDS BY DESIGN. Triphenyltin chloride's InChI
is `3C6H5.ClH.Sn` -- five components before any tool touches it. INPUT
SELECTION took InChI-first unconditionally, inherited from Run_Marvin_pKas.py
where it is correct, so Marvin was handed an already-shattered molecule and
protonated the pieces.

487 compounds have an InChI more fragmented than their SMILES -- cobalamins,
Ni/Fe/Mg porphyrins, molybdenum cofactors, organotins, iron-sulfur clusters.
461 shipped fragmented.

THE RULE is now InChI-first UNLESS the InChI is the more fragmented of the two.
Comparing fragment counts rather than screening for metals keeps it general.
493 compounds take the SMILES route on that test.

Two subtleties, both of which cost a regeneration to find:

  The InChI is DEMOTED, not reordered. Sharing one ladder lets a
  writable-but-wrong rung beat a correct one: on KEGG C18384 the SMILES yields
  the right dative-bonded magnesium propionate, which Marvin's SMILES writer
  refuses, while the disconnected InChI yields a writable three-fragment answer
  and wins.

  The InChI ROW is written from the InChI-derived molecule. Each column carries
  what its own representation can express, which is what 23.4 did: CPD-18407
  ships a connected 8-iron cluster in SMILE and disconnected
  `InChI=1S/C.8Fe.6HS.3S/...` in InChI. This is also load bearing for
  stability. Asking Marvin to write an InChI for a CONNECTED metal cluster
  aborts the JVM with `free(): double free detected in tcache 2` from
  InChINativeGenerateInChICall -- a native fault no handler can catch, and one
  the previous regeneration died on.

RESULT, measured against 23.4 rather than against the previous run:

    fragmented vs 23.4    461 -> 16 -> 8
    R rows                8,727 vs 8,704, 0 invariant violations
    rows                  141,715, exactly 23.4's
    ChEBI/KEGG/Rhea       0 compounds more fragmented than 23.4

The remaining 8 are not input selection. Marvin 26.1 breaks metal coordination
bonds during protonation itself, from either representation: ferrocene splits
into iron and two cyclopentadienyls, C12862 sheds both ammines. Confirmed by
protonating the connected SMILES directly. Recorded, not worked around.

WHY THIS SHIPPED. ModelSEED#296 claimed "0 structures changed" -- true only against
ModelSEED#295, and both were wrong. Comparing a fix to a broken baseline instead of to
23.4 is the check that would have caught it, and is now the first table in the
report.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants